In seismic data, certain geologic features that may not be easy to identify are known to facilitate the discovery of hydrocarbons. These features may be characterized by sub-features or sets of sub-features that are distributed across multiple seismic volumes. For example, seismic reflection data can be used to reconstruct the earth's subsurface structure. Some of these structures may be suitable for trapping fluids such as oil/water or gas. The source to receiver separation-related changes in the amplitude of the reflections from a structure may indicate the type of accumulation (oil/water/gas) in the structure. Currently industry practice is to calculate multi-dimensional (e.g., point, surface, volume, temporal) seismic attributes (i.e., measurements based on seismic data) by using multiple techniques and representations of the data. These attributes aid humans or computer algorithms to interpret and/or identify geological/geophysical features of the data. A severe limitation is that critical attribute interdependencies are not captured during the processes of conventional attribute generation. This limitation arises from the use of different data representations and calculation techniques to generate the attributes; e.g. some attributes are obtained using single trace calculations, whereas others are obtained by operating on multiple traces. (A trace is a record of seismic signal amplitude detected by a single receiver or array of receivers as a function of time.)
Feature Detection Using Current Attribute Extraction Methods
This section reviews conventional approaches to extract attributes. Extraction of attributes from seismic data is a mature area. Local amplitude, phase, frequency, dip, and discontinuity are some of the important, well-studied seismic attributes. Such attributes are then used to facilitate manual or automatic recognition of desired geologic features in seismic data.
The methods used to generate attributes can differ significantly. For example, local frequency, amplitude, and phase attributes are extracted by analyzing each trace separately, whereas dip and discontinuity attributes are extracted by analyzing multiple traces. Due to differences in the attribute generation methods employed, conventional techniques suffer from fundamental limitations including:                1. Conventionally generated attribute volumes can contain a lot of redundant information. For example, peak amplitude and mean-squared amplitude attribute volumes are known to be highly correlated (Barnes, 2006). Such redundant “information” generation not only wastes computational and storage resources, but more importantly overwhelms the end user, thereby resulting in degraded efficiency and missed prospecting opportunities.        2. Each conventional attribute generation method has a different sensitivity to noises in the input datasets. Consequently, the relative quality of generated attribute volumes can differ significantly. Such uncertainty can overwhelm the end user.        
The following sub-sections describe conventional methods to extract routinely used attributes such as local frequency, dip, local amplitude, phase, and discontinuity attributes. In contrast, it will be illustrated later how these attributes can be extracted from the data's curvelet representation by the present inventive method.
Local amplitude and phase attribute: Seismic data typically displays continuity and smooth amplitude and phase content variation across multiple traces. Local amplitude attribute describes how the magnitude of a signal varies across seismic data. Local phase attribute describes how the phase of a signal varies across seismic data. Local amplitude and phase are also referred to as instantaneous amplitude and phase respectively.
The conventional approach computes the analytic (complex-valued) signal separately for each trace in the data set. Such an approach is commonly called complex trace analysis. The magnitude and phase of the computed analytic signal provides the trace's local amplitude and phase attribute respectively. The analytic signal is computed using the one dimensional (1-D) Fourier transform or Hilbert transform.
A key limitation of all conventional local amplitude and phase extraction methods stems from analyzing each trace separately. This approach can create significant additional noise because the individual decompositions are not required to honor the typically smoothly varying amplitude and phase content of seismic reflectors across multiple traces.
Local frequency attribute: Local frequency attribute describes how the local frequency content of a signal varies across seismic data. Local frequency attribute extraction is also referred to as time-frequency representation (“TFRs”), spectral decomposition, or instantaneous frequency. A variety of TFRs such as the short-time Fourier transforms (“STFT”), Wigner-Ville Distribution (“WVD”), adaptive TFRs, complex trace-based methods and Instantaneous Spectral Analysis (“ISA”) have been proposed in literature.
Similar to the conventional local amplitude and phase attributes, conventional local frequency attributes are noisy because the methods analyze each trace separately. An additional limitation of such methods is that the frequency of the seismic data is measured along the direction of a trace. Such a frequency estimate is a useful attribute for flat reflectors. For dipping reflectors, a more useful and physically meaningful attribute is the frequency that is measured in the normal direction to the dip of the reflector because this value more accurately measures the contrast in rock property distribution. To measure the reflector frequency in the normal direction to the dip, it is essential to jointly analyze multiple adjacent traces. Consequently, conventional methods are incapable of accurately measuring the desired local frequency of dipping reflectors (a problem that is solved by the present inventive method).
Dip attribute: This attribute describes the geometric orientation of reflectors in a two dimensional seismic image. Given a reflector in a seismic image, the dip of the reflector at any location is the orientation of the tangent to the reflector with respect to the horizontal. In higher-dimensional seismic volumes, several parameters are necessary to fully describe a reflector's geometric orientation; for example, in three dimensions, the two parameters, commonly referred to as dip and azimuth, are needed to describe a reflector's geometric orientation. For convenience, but without loss of any generality, this document will collectively refer to attributes describing a reflector's geometric orientation as dip.
One class of traditional approaches applies complex trace analysis on individual traces and then calculates the dip. Such approaches can estimate a single average dip at each location. Clearly, such an average dip is not useful when the underlying data contain several events with conflicting dips. Further, such an analysis does not provide the strength of the event containing such a dip. Another class of approaches scans a discrete number of dips at each location in the seismic volume and measures the strength of each dip by computing the local correlation of the data with a plane that has the specified dip. Marfurt (“Robust estimates of 3D reflector dip and azimuth,” Geophysics 71, 29-40 (2006)) further improves upon such dip estimates by performing the dip search over multiple shifted windows. Such approaches can only provide a dip strength that is essentially averaged over all frequencies. Clearly, such average dip strengths are not useful when the underlying data contain events with different frequency content and conflicting dips.
Discontinuity attribute: Discontinuity attribute is closely related to dip attribute because it quantifies rapid changes along the dip direction. Discontinuity attribute is typically used to emphasize fault features or stratigraphic boundaries in seismic data. This attribute is often used interchangeably with the coherence attribute. All conventional discontinuity attributes are based on computing local cross-correlations, as described in U.S. Pat. No. 5,563,949. Such methods are also computationally expensive and can be sensitive to the noises in seismic data.
Background of Enabling Tools
First, this section provides background information about a well known mathematical tool used in the present invention called the curvelet transform. Second, this section reviews known applications of curvelets in the area of geology and geophysics (including one known instance of a curvelet-based attribute extraction), and points out some of their limitations.
Curvelets: A Multi-Resolution Directional Transform
Curvelet transforms are recently developed mathematical tools that represent an m-dimensional signal (m>1) signal using a linear, weighted combination of special elementary functions.
                              signal          ⁡                      (                                          x                1                            ,                              x                2                            ,              …              ⁢                                                          ,                              x                m                                      )                          =                              ∑            i                    ⁢                                    weights              i                        ×                                          elementary_                ⁢                function                ⁢                s                            i                        ⁢                          (                                                x                  1                                ,                                  x                  2                                ,                …                ,                                  x                  m                                            )                                                          (        1        )                                          signal          ⁡                      (                                          x                1                            ,              …              ,                              x                m                                      )                          =                  real          (                                    ∑              i                        ⁢                                          weights                i                            ×              curvelet_elementary              ⁢                                                _                  ⁢                  fns                                i                            ⁢                              (                                                      x                    1                                    ,                  …                  ⁢                                                                          ,                                      x                    m                                                  )                                              )                                    (                  1          ⁢          a                )            For example, FIG. 1 illustrates a synthetic image, which comprises two intersecting reflectors and its representation as a sum of weighted curvelet elementary functions. For illustration purposes, only four curvelet elementary functions are shown. Alternatively, curvelet transforms can also express an m-dimensional signal as
                              signal          ⁡                      (                                          x                1                            ,              …              ,                              x                m                                      )                          =                  real          ⁡                      (                                                                                                      ∑                      i                                        ⁢                                                                                            complex_                          ⁢                          weights                                                i                                            ×                                                                                                                                        complex_elementary                    ⁢                                          _fns                      i                                        ⁢                                          (                                                                        x                          1                                                ,                        …                        ⁢                                                                                                  ,                                                  x                          m                                                                    )                                                                                            )                                              (        2        )            Due to the properties of the elementary basis functions (similar to real and imaginary parts of an analytic signal), the complex valued representation in Equation (1) can be transformed into Equation (2) and vice versa with simple manipulations on the equation weights. For example, the implementation in the Curvelab Software Package, available for download in Matlab and C++ at www.curvelet.org (2005), provides an option to compute the complex-valued curvelet transform, which is equivalent to the representation in Equation (2).
Properties of curvelet elementary functions: Curvelet transforms provide desirable representations for seismic data due to multi-resolution and directional properties of the elementary functions. As illustrated in FIG. 2, the elementary functions are highly oscillatory in one direction (referred to as the oscillatory direction henceforth), with the function's cross-section 21 in the oscillatory direction resembling a windowed sinusoid. Since the oscillations in different elementary functions can occupy different frequency bands (which is a resolution measure), the transform is said to possess multi-resolution properties. The transform is said to possess directional properties because the elementary functions are not isotropic. They vary more slowly (i.e., contain slower frequency oscillations) and smoothly in the directions orthogonal to the oscillatory direction; the function's cross section 22 in the orthogonal direction resembles a Gaussian window. Each elementary function is spatially localized because its amplitude rapidly decays to zero outside a certain region. The region that essentially localizes an elementary function resembles an oriented needle in two dimensions (2-D) and a disc in three dimensions (3-D). The location, oscillation frequency, and oscillatory direction of each elementary function can be decoded from its index i in Equation (1) or (2) because the functions are ordered. Due to the listed properties of the elementary functions, the curvelet transform belongs to the class of multi-resolution directional transforms (see, for example, U.S. Patent Application Publication US 2007/0038691).
Obtaining curvelet representations (Equations (1) or (2)): A given dataset's curvelet representation is not unique because the number of curvelet elementary functions exceeds the number of points in seismic data. Various curvelet representations differ in their choices of weights in Equation (2).
The most common curvelet representation is obtained by taking the forward curvelet transform of a given dataset. The forward curvelet transform can be implemented by employing a variety of steps that will be familiar to one versed in the multi-resolution directional transforms area (see, for example, Candes, et al., “Fast discrete curvelet transforms,” SIAM Multiscale Model. Simul. 5, no. 3, 861-899 (2006, published online in July 2005)). Some key steps comprise taking the multidimensional Fourier transform, windowing the Fourier coefficients into annuli (termed subbands or scales sometimes) with dyadic widths, further subdividing the annuli into angular wedges, and finally taking inverse Fourier transforms (see the Candes, et al. article referenced in the preceding sentence for details).
An alternate curvelet representation that employs a sparse set of weights in Equation (2) can be obtained by employing minor variations of the iterative method described in Daubechies et al., “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” Communications on pure and applied mathematics 57, 1413-1457 (2004):                a. Compute an initial estimate of the weights by employing the forward curvelet transform.        b. Process the weights by shrinking (see hard thresholding and soft thresholding, for example, in Mallat, A Wavelet Tour of Signal Processing, Academic Press, 441-446 (1998)) the weights, thereby leaving the large weights relatively untouched.        c. Compute the error between the original signal and the signal represented by the processed weights.        d. Compute the weights for the error signal (similar to step (a)).        e. Add the error signal weights to the processed weights in step (b) to create an updated set of processed weights.        f. Go to step (b) and iterate.        
Curvelet representations with a sparse set of weights can also be computed by employing iterative techniques such as Matching Pursuit and Orthogonal Matching Pursuit (see Tropp and Gilbert, “Signal recovery from random measurements via Orthogonal Matching Pursuit.” Submitted for publication, April 2005. Revised, November 2006. Download: (http://www-personal.umich.edu/˜jtropp/papers/TG06-Signal-Recovery.pdf); Basis pursuit (see the previously mentioned Mallat, A Wavelet Tour of Signal Processing, pg. 409); LASSO (see Tibshirani, “Regression shrinkage and selection via the LASSO,” J. Royal. Statist. Soc B. 58(1), 267-288 (1996)); LARS (see Efron et al., “Least Angle Regression”, Ann. Statist. 32(2), 407-499 (2004)); and expectation-maximization, Bayesian estimation algorithms, belief propagation or similar techniques referenced in these publications.
Alternatives to curvelet transforms: The discussion above has dealt with the curvelet transform. However, several other multi-resolution directional transforms such as the 2-D and 3-D complex wavelet transform (Selesnick et al., “The Dual-Tree Complex Wavelet Transform,” IEEE Signal Processing Magazine 22, 123-151 (2005)); the 2-D contourlet transform (Do and Vetterli, “The contourlet transform: an efficient directional multiresolution image representation,” IEEE Transactions Image on Processing, 14, 2091-2106 (2005); Nguyen and Oraintara, “A Shift-Invariant Multiscale Multidirection Image Decomposition,” 2006 IEEE International conference on acoustics, speech, and signal processing. Volume 2); and the 3-D surfacelet transform (Lu and Do, “Multidimensional Directional Filter Banks and Surfacelets,” IEEE Transactions on Image Processing 16, 918-931 (2007)) provide data representations that share the properties and strengths of curvelet representations. Steerable pyramids (Freeman W. and Adelson R., 1991, “The design and use of steerable filter,” IEEE Pattern Analysis and Machine Intelligence 13, 891-906.) are another example of transforms that possess multi-resolution and directional properties. These alternate multi-resolution directional transforms can be substituted for curvelet transforms in some applications.
Currently available implementations for curvelets (for example, Curvelab Software Product) make fixed choices for the windows, the dyadic widths of the annuli, and number of wedges based on theoretical considerations. A person skilled in the area of multi-resolution directional transforms would recognize that a curvelet-like multi-resolution directional transformation can also be implemented by varying the windows employed (e.g., use Hamming windows instead of Meyer windows), the annuli widths (e.g., not necessarily dyadic), and the number of wedges. For example, the previously mentioned contourlets in 2-D and surfacelets in 3-D choose the windows, annuli, and wedges differently from the curvelet transform (Curvelab).
Further, current implementations such as the Curvelab product employ the multidimensional Fourier transform to compute the forward curvelet transform. A person skilled in the area of multi-resolution directional transforms would recognize that curvelet-like multi-resolution directional transformations could also be implemented via equivalent computational operations in the data domain. For example, the previously mentioned complex wavelet transform, the contourlets in 2-D, and the surfacelets in 3-D can all be implemented by performing computations in the data domain.
Applications of Curvelets and Related Transforms
Geophysical data processing: Curvelets have found several successful applications in enhancing the resolution of a geophysical dataset. Recently, Herrmann and Verschuur (“Curvelet-domain multiple elimination with sparseness constraints,” 74th SEG Annual Meeting, Expanded Abstracts (2004)) and Yarham et al., “Curvelet-based ground roll removal,” 76th SEG Annual Meeting, Expanded Abstracts, SPNA 1.7 (2006)) attenuated undesirable water bottom multiples and ground roll noise, respectively, from seismic data by appropriately attenuating the data's curvelet transform coefficients. Moghaddam and Herrmann (“Migration preconditioning with curvelets,” 74th SEG Annual Meeting, Expanded Abstracts, 2204-2207, (2004)) and de Hoop and Douma (“On common-offset prestack time migration with curvelets,” 75th SEG Annual Meeting, Expanded Abstracts, SPMI P1-3 (2005)) employ curvelet transforms to improve the resolution of migrated seismic data. Hennenfent and Herrmann (“Sparseness-constrained data continuation with frames: Applications to missing traces and aliased signals in 2/3-D,” 75th SEG Annual Meeting, Expanded Abstracts (2005)) use curvelets to recover higher resolution, interpolated geophysical data. These curvelet-based data enhancement methods are unrelated to the technical problem addressed herein because they apply curvelet transform to problems that are fundamentally different from the identification of geologic patterns from seismic data.
Curvelet-based seismic texture classification: Two papers (Monsen and Odegard, “Seismic texture analysis by curve modeling with applications to facies analysis,” SEG Expanded Abstracts 21, 617 (2002); and Monsen et al., “Seismic Texture Classification by Hidden Markov Tree Modeling of the Complex Wavelet Transform,” Norwegian Signal Processing Symposium, Trondheim, Norway (October 2001)) describe the limited use of complex wavelets and contourlets respectively to calculate two specific textural attributes (parallel and chaotic) of a two-dimensional seismic image. The textural attributes are then used to classify a two-dimensional image. The classifications presented in both papers are based on the statistics of the seismic data's complex wavelet and contour representation weights. Monsen et al. (2001) use a hidden Markov model on the data's complex wavelet coefficients to perform the classification. Monsen and Odegard (2002) base their classification on a simpler statistical model, which assumes that the data's contourlet representation weights in different sub-bands and angular wedges are independent. The authors speculate that the technique could be extended to 3-dimensional problems.
PCT International Publication No. WO 2006/064239 A1 discloses a complex wavelet-based method to identify features within a dataset. The method uses a complex-valued multiscale transform such as the complex wavelet transform. This work applies a complex-valued multiscale transform to a target dataset and a candidate dataset, and then uses the phase differences between coefficients of the two datasets to generate a measure of match. This measure is subsequently used to accept or reject a hypothesis that the target is present. No attributes or their interdependencies are extracted or retained from the complex-valued multiscale representation of the data.
In a patent application publication (US 2007/0223788 (pub. date Sep. 27, 2007); CA 2571094 (pub. date Jun. 13, 2007)), Pinnegar et al. disclose a method to compute frequency-dependent dip or amplitude maps by processing multidimensional data using space-frequency or time-space-frequency transforms. The maps are constructed from dominant features that correspond to the largest value of the transform coefficients. These attributes are created to assist in the interpretation of geologic features from seismic data. The authors suggest that the curvelet transform (a space-frequency transform) can be used to create the attribute maps.
In another patent application publication (US 2007/0043458 (published on Feb. 22, 2007); WO 2007/006145 (published on Jan. 18, 2007)), Pinnegar discloses a method for processing time-varying three component signals to determine polarization dependent attribute maps. The attribute maps are determined by transforming each of the three components to a time-frequency domain. The author focuses on using the S-transform to construct the attribute maps, and suggests that the curvelet transform can be used instead of the S-transform.
In a Ph.D. thesis (“Seismic applications of complex wavelet transforms,” Delft University of Technology, 2002, ISBN number 90-9015810-3) and another publication (“Non-redundant directionally selective complex wavelets,” 2000 International conference on image processing, pages 379-382), van Spaendonck describes the use of a complex wavelet transform to compute dip and frequency (scale) seismic volume attributes. In another paper (“Directional scale analysis for seismic interpretation”, SEG Expanded Abstracts 1844-1847 (1999)), van Spaendonck describes the application of steerable pyramids to compute dip and frequency (scale) seismic volume attributes.
What is needed is an expansive use of a more powerful transform such as the curvelet transform to develop new seismic data attributes that the user may decide upon (not limited to specific attributes such as the chaotic and parallel textural attributes) such as (without limitation) local frequency, dip, local amplitude, phase, and discontinuity, and preserve interdependencies between multiple attributes. The present invention satisfies this need.